clear all
path(path,'../TOOLS_01');

[parametros_malha] = le_malha('Malha_problema_inverso_5_nos.msh');

n_eletrodos=parametros_malha.n_eletrodos;
n_elementos=parametros_malha.n_elementos;
n_nos=parametros_malha.n_nos;
n_objetos=parametros_malha.n_objetos;
coords=parametros_malha.coords;
topologia=parametros_malha.topologia;
topologia_eletrodo=parametros_malha.topologia_eletrodo;
idx_malha=parametros_malha.idx_malha;
idx_objeto=parametros_malha.idx_objeto;

%parametros gerais
altura=0.05; % altura da malha
rho_eletrodo=1/50; % resistividade do eletrodo
corriente=1e-3; % corrente aplicada


%%%% valor inicial

rho=2.0;
rho_objeto=rho;

rho_malha=zeros(1,n_elementos+n_eletrodos);
rho_malha(idx_malha)=rho;
rho_malha(idx_objeto)=rho_objeto;
rho_malha(n_elementos+1:end)=rho_eletrodo;
alpha=1e-9;
%%%% valor inicial


[ centroides ] = calc_centroides( coords, n_elementos,topologia);
max_x=max(centroides(:,1));
min_x=min(centroides(:,1));
max_y=max(centroides(:,2));
min_y=min(centroides(:,2));
res=64;
[XI,YI]=meshgrid(linspace(min_x,max_x,res),linspace(min_y,max_y,res));

%%%% matriz regurarizacao
%%% calcula o raio maximo
R_MAX=-1.0e9;
for i=1:n_eletrodos
    r=norm(coords(topologia_eletrodo(4),:));
    if (r > R_MAX)
      R_MAX=r;
    end
end
%%% calcula o raio maximo
norm_value_FtF=30.0; % desvio padrao do filtro gaussiano ATENCAO: valor em porcentagem do raio maximo 
norm_value_FtF=(norm_value_FtF/100.0)*R_MAX;
[FtF]=calc_FtF_gaus(centroides,n_elementos,alpha,norm_value_FtF);




load lixo.mat

n_iteracoes=30;
passo_nr=1;

for i=1:n_iteracoes
    [Y_global,Y_local] = monta_rigidez_tri(n_nos,n_elementos,coords,topologia,rho_malha,n_eletrodos,topologia_eletrodo,altura);
    [ H,V_calculado ] = monta_jacobiano_2D(Y_global,Y_local,rho_malha,corriente,n_elementos,n_eletrodos,n_nos,coords,topologia_eletrodo);
    delta_rho=(H'*H+FtF)\(H'*(V_objeto-V_calculado)'-FtF*rho_malha(1:n_elementos)');
    rho_novo=rho_malha(1:n_elementos)+0.5*delta_rho';
    rho_novo(find(rho_novo<=0))=1.0;
    %%%%% plota imagem
     F = TriScatteredInterp(centroides(:,1),centroides(:,2),rho_novo');
     QZI=F(XI,YI);
     figure(2),surf(XI,YI,QZI),title('\rho'),view(0,90),colorbar
     
     F1 = TriScatteredInterp(centroides(:,1),centroides(:,2),delta_rho);
     QZI1=F1(XI,YI);
     figure(3),surf(XI,YI,QZI1),title('Delta \rho'),view(0,90),colorbar
     
%      pause
    %%%%% plota imagem
    rho_malha(1:n_elementos)=rho_novo';
    
    erro(i)=norm((V_calculado-V_objeto));
    if (i>4)
        if erro(i)>erro(i-3)
            passo_nr=passo_nr;
        end
    end
            
    figure(1),semilogy(erro,'.-'),grid,title(sprintf('passo NR = %g',passo_nr))
end


F = TriScatteredInterp(centroides(:,1),centroides(:,2),rho_novo');
QZI=F(XI,YI);
figure(20),surf(XI,YI,QZI),colorbar


figure(1),semilogy(erro,'.-'),grid



